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In this paper we consider the possibility to use numerical simulations for a computer assisted 
analysis of integrability of dynamical systems. We formulate a rather general method of recovering 
the obstruction to integrability for the systems with a small number of degrees of freedom. We 
generalize this method using the results of KAM theory and stochastic approaches to the families 
of parameter depending systems. This permits the localization of possible integrability regions in 
the parameter space. We give some examples of application of this approach to dynamical systems 
having mechanical origin. 
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I. INTRODUCTION 

This paper is the first in a series of works devoted to description of the possibilities of application of numerical 
methods to analysis of integrability of dynamical systems. The problem of integrability has been studied since the 
middle of the XIX century, when the question of primary interest was to be able to integrate the system of differential 
equations by inversion of functions and quadratures, that is give a more or less explicite solution. Nowadays with the 
development of methods relating symplectic geometry and dynamics the notion is extended, namely one studies the 
existence of an appropriate number of conserved quantities (first integrals, invariant measure, etc), possessing some 
properties. 

In this paper we will consider integrability in the Liouville-Arnold sense, namely for a hamiltonian system with 
n degrees of freedom it is the existence of n independent first integrals in involution ([!]). In particular we will 
pay attention to the independence condition from the topological point of view. We propose a constructive method 
of computer assisted proof of non-integrability of systems with small dimensional configuration space based on the 
analysis of the topology of the phase space. We generalize it using the results of the Kolmogorov- Arnold- Moser theory 
to the systems with parameters. We also provide examples of application of this method to some dynamical systems 
having mechanical origin. 



II. METHOD OF SECTIONS 



Consider an autonomous hamiltonian system with a two-dimensional configuration space Q, its phase space T*Q is 
of dimension 4. Since the right hand sides of the equations of motion do not depend explicitely on time, any trajectory 
of this systems belongs to a constant hamiltonian (energy) level hypersurface having the dimension 3. For complete 
integrability of the system with two degrees of freedom the existence of another first integral independent with energy 
integral is needed - in this case any trajectory would belong to the 2-dimension manifold defined by the intersection 
of the level surfaces of these first integrals. 
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A. Idea of the method 



This simple topological consideration alone does not give a method of analysis of integrability, but together with a 
good visualization algorithm it permits to give an answer for a rather large class of systems. Effective visualization of 
the dynamics in a 4-dimensional phase space is not an easy task, therefore we consider the intersection of a trajectory 
of the system with 2-dimensional planes in it. Let us discuss the possible results of this intersection. 

The dimension of intersection of two generic manifolds of dimensions 71 1 and 712 in an A^-dimensional space is given 
(0) by the equation 

dim — Til + n2 — N. (1) 

In our case the manifold swept by the trajectory intersects a 2-dimensional plane (712 — 2) in the phase space of 
dimension TV = 4. Energy conservation guarantees that the dimension (ni) of this manifold is at most 3. Thus we can 
observe two possible cases depending on the existence of additional first integral: either rii — 2 and the intersection 
is 0-dimensional (finite set of points) , or ni = 3 and the intersection is of dimension 1 (finite set of curves) . It is clear 
that the presence of curves in the intersection is an obstruction to complete integrability. Abscence of curves however 
does not directly mean integrability, since it can be the consequence of various reasons - we will discuss them more 
precisely in the section [Till in the context of generalization of this method. 

Let us be more concrete on this idea about intersection and prove in the particular case the formula ([T]). Choose 
the canonical coordinates qi,pi,q2,P2 of the phase space T*Q (assume for simplicity that it can be done globally). 
Then an arbitrary two-dimensional plane (or, better to say two-dimensional linear subspace) is given by the system 
of equations 

aigi + hpi + CiQ'2 + dip2 = ei, 

a2qi + b2Pi + C2q2 + d2P2 ^ £2- (2) 

for some constants a^, 6^, q, di, e^, z = 1, 2. Since the system is autonomous any trajectory belongs to the energy level 

h=H = h. (3) 

The conditions ([2l[3]) on 4 coordinates in the phase space define a 1-dimensional manifold. If there exists another first 
integral 

I2 = const (4) 

the system ([21 El H]), if it is compatible, admits a finite set of solutions. It is now easy to see that it is important to 
pay attention to the independence of the conditions ([2]) from Q and This is not difficult to guarantee for the 
energy integral, since it is non-linear and usually known; if the other first integral exists and the equation Q turnes 
out to be dependent with ^ it means that the system admits some degeneracy, that will be seen in the numerical 
simulation. To avoid "false detection" by this method it is enough to consider the intersection with a couple of 
mutually independent planes, that is done in practice. 



B. Triple pendulum 

As one of the main examples of application of the method in this work we will consider the flat motion of the 
pendulum-type systems. A multiple pendulum is the system of point masses, connected by weightless inextensible 
rods, the first of the points being fixed. For a triple pendulum these conditions read 

(r, -r,_i)2 = Z, i = l,2,3, 

where corresponds to the i-th mass, the fixed point is Tq. One can also consider a more general case of constraints 
given by arbitrary polynomials of degree 2, which however does not always correspond to a physical configuration. 

The choice of this system is motivated by several factors. Let us note that the case of a double pendulum is 
well studied. In particular its free motion on the plane is a classical example of a completely integrable system 
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(see for example [3j). In the notation of the previous section taking the angular momentum for I2 one obtains two 
independent conditions ([31 H]) , having a 0-dimensional intersection with ^ . But already in the presence of gravity 
the non-integrability can be shown and some chaotic behaviour has been observed, that is the trajectory is rather 
dense on the energy level-surface, that can be seen on the sections. The problem of control has also been studied for 
this system (Q)) among interesting results one can mention that a double pendulum can be stabilized in the upright 
position by controlling only one degree of freedom. 

Meanwhile the dynamics of the triple pendulum is almost not studied, although it is interesting for applications. We 
have shown rather directly (without going much into details about topology, but writing explicitely the coordinates) 
that the behaviour of the system is rather irregular and it is not integrable |173 • This results can be useful to study 
the dynamics of a multiple lattice, being one of the low-temperature models of the polymer molecule. The idea of 
describing the interaction of subunits in molecular dynamics by lagrangian mechanics extends the harmonic model 
used in [5!] . Actually the study of thermodynamic properties of the pendulum systems in the framework of the model 
from 5] explains the author's interest to their integrability. 

We will not describe this example in full details since we have already done it partially in ^ and in Q. Let us 
just note that in the original setting the system obviously has 3 degrees of freedom, but since it admits the angular 
momentum first integral it can be reduced by the Routh transform (Q) to the one having the 4-dimensional phase 
space, that is precisely the range of application of the method of sections. Let us also note, that this reduction is 
profitable to understand the topology of the phase space, but to perform numerical simulations it is more convenient 
to use the lagrange multipliers (Q) and then compute the reduced coordinates from the cartesian ones. 

A typical result of the numerical simulations is shown on fig. [T] The coordinates in the reduced phase space are the 
angles /3i,/32 between the segments of the pendulum and their derivatives. In the hamiltonian formalism one should 
actually use the momenta instead of velocities, but we use the natural duality of T*Q and TQ, so it gives an equivalent 
picture from the topological point of view. The intersecting planes are chosen to be parallel to the coordinate ones. 
The presence of curves in the intersection clearly shows non-integrability. 




FIG. 1: The intersection of the trajectory of the triple pendulum with the planes; 

a). (/3i = l,/32 = 1), b). (/3i,/32), c). (/Ji,/3i), d). {132, P2). 
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C. Satellite dynamics 



Let us consider another mechanical example for which the method is applicable: the motion of the dynamically 
symmetric satellite on the circular orbit. In the orbital coordinate system, the axes of which are directed along the 
radius of the orbit of the center of masses of the satellite, its normal and its binormal, the position of the satellite is 
given by the three Euler angles "0, 0, Lp. Denoting the corresponding momenta by p^^pg,p^, one obtains the following 
hamiltonian function: 

^ = o ^"^2 n + ^ ~ Pipctg{e) cos(V') - aPp^^^^ - pe sin(V') + 
2 sm a I sm i) 

+a/3^ + — ^ + - a - 1 cos2 0, 5 
sm 2 sm 9 2 

where a = C/A; A,B,C are the principal moments of inertia {A = B). The coordinate (f is cyclic, thus we can 
fix the corresponding momenta p^p = a(3 = const, where /3 denotes the ratio of the orbital angular velocity and the 
projection of the angular velocity of the satellity to its symmetry axis. Let us consider the values of a = 4/3,/? = 0, 
for which the stability of equilibrium solutions is studied in [l^l . The reduced system has two degrees of freedom and 
is described by the hamiltonian: 

2 2 1 

2 sm 2 2 

That is we can again apply the method of sections. The intersections that one obtains (fig. ^ show its non-integrability. 



Let us however note, that to obtain the sections with the curves for this problem we had to analyse a lot of different 
initial conditions. This indicates that the system can be locally integrable, that is possess an additional first integral 
for a certain subset of initial data. 



D. Details of the method of sections 



Let us now discuss the possibilities of application of the method of sections to other systems, in particular let us 
comment on the class of systems for which it can result in the precise conclusion on integrability. As we have noted 
in the very beginning the natural restriction comes from the dimension of the phase space, or more precisely from the 
possibility to reduce it. This however covers a large class of systems for which the integrability is an open question, 
such as the dynamics of the triple lattice, the motion of a mass point in the symmetric potential in a 3-dimensional 
space, dynamics of geodesies on the curved surfaces etc. 

The only difficulty is the choice of convenient coordinates for constructing the sections. But this problem can be 
naturally solved when the system admits the angle-parametrization like in the described examples. Note that such 
systems often arise in the applications. According to the theorem by V.V.Kozlov ([11]) it is interesting to consider the 
systems with the genus of the configuration space at most 1, since if it is not the case one knows that the system is 
not analytically integrable; but it means that basically any choice of coordinates has the structure of angles. We can 
however study more general surfaces if we are interested in a larger class of first integrals. Then it is good to make 
sure that the motion takes place in a bounded domain of the phase space, to accumulate a sufficiently representative 
intersection. This condition is also rather natural at least for the systems with compact configuration space. 

Let us also note a technical difficulty of the method. First, to conclude non-integrability one needs to compute 
a rather long trajectory, since the intersection of the trajectory with a surface of small dimension is a rear event in 
contrast, say, to Poincare sections where every loop gives a point in the plane. But this reduction of the dimension is 
done intensionally, as it permits more effectively classify different cases: it is much easier to distinguish points from 
curves than curves from thin domains that can occur in the Poincare sections. Our method also permits not to pay 
much attention to transversality of the trajectory to the planes. All this only results in the need to apply reliable 
algorithms of numerical integration. 
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FIG. 2: The intersection of the trajectory with the planes, parallel to the coordinate ones, containing the point t/) = 0.28, Q = 
0.82, = 0.15, pe = 0.37. 

Another thing that we haven't yet discussed in details is the interpretation of the "empty" (O-dimensional) inter- 
sections. The existence of such sections does not permit to make any direct conclusion since it only means that a 
given trajectory is by chance more regular. But in the next section we will explain the origin of such sections for 
non-integrable systems and discuss the possibility to use them to qualitatively describe the behaviour of the systems. 



III. GENERALIZATION OF THE METHOD OF SECTIONS VIA THE KAM THEORY 



The method of sections as we have seen permits to prove non-integrability for a given dynamical system. This 
problem however is rather special, more often one wants to find the relation between the parameters when integrability 
is possible. In this section we will consider an example of a two parameter system showing how one can extend the 
method of section to this problem. Despite the discrete nature of the method it permits (via some extra mathematical 
considerations) also to draw conclusions for a continuous range of the parameters. 



A. Pendulum- type systems 

Consider free flat motion of a system obtained from the triple pendulum by moving the fixing points of the next 
segment along the previous one (fig. [3]). We shall call such objects pendulum- type system. In particular they describe 
the motion of a physical pendulum, i.e. rigid bodies fixed between themselves. 
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FIG. 3: A pendulum type system parametrized by the angles. The distances: El = eih, E2 = £2/2 



The configuration of this system is described by two dimensional vectors ri and velocities Vi {i = 1, 2, 3) with the 
following constraints: 

^1 = (ri - ro f -ij^O, 
(P2 = (r2 -eiri f -l^^O, 
f3 = (rs - (eiri + £2(r2 - £iri)))^ - ^3 = 0, (6) 

the point rg is fixed. The parameters Ei characterize linearly the fixing point of the {i + l)-st segment to the z-th one: 
Si = 1 corresponds to fixing at the endpoint (like in the multiple pendulum), = - to the fixed point (like two 
non-interecting pendula with the same fixed point). 

Integrability. Note that for ei — 62 — I v/e obtain the non-integrable case of the triple pendulum studied before. 
For £1 = and any 62 the system decouples to non-interacting simple pendulum and a double pendulum-type system, 
and becomes obviously integrable. And in the particular case when both £1 = £2 = 0, it decouples to three non- 
interacting simple pendula. For other values of Si we will be again interested in the existence of a sufficient number 
of independent first integrals. 



B. Topology of the phase space 



For all values of ei, £2 the system is invariant under rotation around the fixed point Vq, thus it possesses the Nother 
integral of angular momentum. Since the position is still characterized by the three angles (fig. [3]), this integral permits 
to perform Routh reduction to the 4-dimensional phase space. To make this reduction explicitely we again study the 
dynamics in terms of the angles between the segments of the pendulum- type system /3j — ai+i ~ cti, i — 1,2 and 
their derivatives. Now for any fixed couple of (£i,£2), one we can apply the method of sections. 



Typical results for the pendulum-type systems are represented on figs. 4(a) 4(d) showing the intersection with the 
planes 02 = 0, /32 = 0), 02 = 0,^2 = 0), 0i = 0,^2 = 0) 02 = iji = 1), 02 = 0, A = 0), 0i = 0, /3i = 0) 

This figures show that the sections contain curves, it means the corresponding systems do not admit an additional 
first integral. The similar figures have been obtained (i.e. initial conditions found - [Tsj ) for all couples £i,£2 G 
[0, 1] X [0, 1] with the step 0.02, except the trivial ones (£1 = 0), that means non-integrability of the corresponding 
systems. 



(a) ei = 0.7, £2 = 0.8 



(b) ei = 0.72, £2 = 0.74 



(c) £i = 0.72, £2 = 0.76 (d) £i = 0.74, £2 = 0.78 



FIG. 4: Sections for various values of ei, 



C. Generalization to continuous parameters 



Note that the sections for close values of parameters and close initial data look rather similar (Fig. |4(b)[[4(c) ). One 



thus is tempted to conjecture that the systems between the two values of parameters are also non-integrable. Alone 
this statement is certainly false, but we can make sense out of it by some extra consideration. Let us make it more 
precise. 

The results of Poincare ([13, show that the systems obtained from integrable one by a perturbation usually fail 
to be integrable. However according to the Kolmogorov-Arnold-Moser theorem (see for example [l|) for the systems 
close to integrable a set of positive measure of invariant tori preserves its topology. In our terms it means, that when 
the initial point even of a non-integrable system lies on such a torus the sections do not contain curves. And when 
we approach an integrable system in the space of parameters the probability to be on such a torus increases. 

This effect is indeed observed in the neighborhood of the integrable case (ei = 0) for arbitrary chosen initial points. 
The figure 5(a) shows the distribution obtained by the Monte-Carlo method ([ij|) of such tori, i.e. the proportion 



of "empty" sections if one starts from a pseudo-random point in the phase space for the given values of parameters 



(ei,£2). The figure 5(b) shows the same distribution only depending on £\. One sees from them, that except the 
segment e\ = 0, in the square (ei,£2) G [0, 1] x [0, 1] there are no pronounced maxima, that shows non-integrability 
of all the systems with parameters different from £\ — 0. That is in the family of triple pendulum-type systems the 
only integrable situation corresponds to decoupling the systems to non-interacting subsystems of smaller dimensions. 



Generalization of the method of sections. It is clear that the same approach is not limited to the above 
example. We can apply it to any family of systems with arbitrary number of parameters the phase space of which 
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(a) Histogram depending on (ei,£2) G (0, 1) X (0, 1). Scale: from — (b) Histogram depending only on ei G (0, 1) 

white, to 16 — black. The maxima are concentrated along the 
segment £i = 

FIG. 5; Distributions of the Liouville tori preserving their topology 

can be uniformcly reduced to the 4-diniensional. One needs to choose a domain D in the parameter space, apply 
the method of sections to a sufficient number of pseudo-random initial points for all values of parameters from a 
rather dense set in D and compute the estimation of the probability to start with an invariant torus its preserving 
the topology. The subdomains of D corresponding to maxima of the obtained distribution are the candidates for 
integrability. 



IV. CONCLUSION 



In this paper we have described the topological properties of the phase space of the small dimensional dynamical 
systems in the form convenient for verification. We have proposed a method of computer assisted proof of non- 
integrability for a given dynamical system and generalized it to the continuous range of parameters. This is a natural 
application of the classical Monte-Carlo method, permitting by a well-developped technique, using the pseudo-random 
quontities construct a qualitative picture of the behaviour of a deterministic system. Together with the results of the 
Kolmogorov-Arnold-Moser theory it permits to answer rather complicated questions on integrability. Such statistical 
approaches are often used in quantum mechanics and in various inverse problems of image recognition. 

As we have mentioned in the beginning of this paper there are other ways to study dynamical integrability using 
numerical methods. For example the meromorphic non-integrability of the examples discussed in the section HIl can 
be proved by constructing the monodromy group ([l^) and application of the results of ([IB])- One of the main 
motivations to develop these methods for us is studying more involved problems like systems with delay naturally 
appearing in biological modeling or relativistic celestial mechanics. 
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